function DecodingProbability_v24()

% Clear workspace
clear all
close all
clc

% Time it
tic

% Transmitter range
T_min=10;
T_step=1;
T_max=70;

% Field cardinality
q=2^8;

% Layer selection probability
g_1=[0.4 0.6]; % 0.5 0.5
g_2=[0.3 0.7]; % 0.5 0.5
g_3=[0.5 0.5]; % 0.5 0.5

% Layer dimensions (these are dejan examples)
k1=15;
k2=25;
K1=k1;
K2=k1+k2;

% Calculate decoding probabilities #set 1
l1_prob_1=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g_1);
l2_prob_1=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g_1);

% Calculate decoding probabilities #set 2
l1_prob_2=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g_2);
l2_prob_2=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g_2);

% Calculate decoding probabilities #set 2
l1_prob_3=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g_3);
l2_prob_3=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g_3);


% Plot nice graph!
hold all

% Linespec options here
% http://www.mathworks.se/help/techdoc/ref/linespec.html

legend_h=plotter(l1_prob_1,l2_prob_1,T_min,T_step,T_max,K1,K2,'r',':','d');
legend_h=plotter(l1_prob_2,l2_prob_2,T_min,T_step,T_max,K1,K2,'b',':','s');
legend_h=plotter(l1_prob_3,l2_prob_3,T_min,T_step,T_max,K1,K2,'black',':','o');

% Legend fix
label_1=strcat('L1, \Gamma=',num2str(g_1(1)),', ',num2str(g_1(2)));
label_2=strcat('L2, \Gamma=',num2str(g_1(1)),', ',num2str(g_1(2)));
label_3=strcat('L1, \Gamma=',num2str(g_2(1)),', ',num2str(g_2(2)));
label_4=strcat('L2, \Gamma=',num2str(g_2(1)),', ',num2str(g_2(2)));
label_5=strcat('L1, \Gamma=',num2str(g_3(1)),', ',num2str(g_3(2)));
label_6=strcat('L2, \Gamma=',num2str(g_3(1)),', ',num2str(g_3(2)));
set(legend_h,'String',[label_1;label_2;label_3; label_4;label_5;label_6],'location','SouthEast')
set(legend_h,'interpreter','tex')

% Save plot
print(gcf,'-depsc2','uep_ew_analytic.eps')

% Time it
toc

end

% Calculate layer 1 probability
function l1_prob = layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

l1_prob=zeros(T_max,1); % Decoding 1. Layer probabilities
sol=zeros(T_max+1,1); % Array for temporary results

for tx =T_min:T_step:T_max % For each number of recv packets
    
    for n = 0:tx % For each permutation of a number of recv packets
        
        % Layer 1 by itself
        val=PM(n,K1,K1,q);
        
        val_test=0;
        
        % Layer 1 and Layer 2 gives rank K2 (This way we also get L1!)
        % Sum prob for all possible ways to achieve rank K2 with given permutation of recv packets
        
        for i=0:K1-1
             tmp1=PM(n,K1,i,q);
            
            %Optimization We only need to calculate tmp2 if tmp1!=0
            if tmp1~=0
                tmp2=PM(tx-n,K2-i,K2-i,q);
                val_test=val_test+tmp1*tmp2;
            end
            
            
        end
        
        % We have counted all the ways L1 can become full rank by itself
        % We have counted all the ways L2 can become full rank (except when
        % L1 is full rank!)
        % Both outcomes are valid for getting L1 and since they are disjoint we can just add them!
        % P(a)+P(b)-P(ab), where P(ab)=0 because they are disjoint!
        
        sol(n+1)=(val+val_test)*binopdf(n,tx,g(1));
        
    end
    
    l1_prob(tx)=sum(sol);
    sol=zeros(T_max,1);
    
    disp(['Layer 1: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Calculate layer 2 probability
function l2_prob = layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

% Decoding 2. Layer probabilities
l2_prob=zeros(T_max,1);

% Array for temporary results
sol=zeros(T_max+1,1);

for tx =T_min:T_step:T_max % For each number of recv packets
        
    for n = 0:tx % For all permutations of recv packets
        
        val=0;
        
        for i=0:K1 % For all possible ways to achieve rank K2 with given permutation of recv packets
            
            tmp1=PM(n,K1,i,q);
            
            % Optimization no need for tmp2 when tmp1=0!
            if tmp1==0
                continue;
            end
            
            tmp2=PM(tx-n,K2-i,K2-i,q);
            val=val+tmp1*tmp2;
        end
        
        sol(n+1)=val*binopdf(n,tx,g(1));
        
    end
    
    l2_prob(tx)=sum(sol);
    sol=zeros(T_max,1);
    
    disp(['Layer 2: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Plotter for a nice graph!
function legend_h = plotter(l1_prob,l2_prob,T_min,T_step,T_max,K1,K2,color,style,symbol)

% Replace 0 with NaN in (l1_prob,l2_prob) for prettier plot
for k=1:length(l1_prob)
    
    if l1_prob(k)==0
        l1_prob(k)=NaN;
    end
end
for k=1:length(l2_prob)
    if l2_prob(k)==0
        l2_prob(k)=NaN;
    end
    
end

% Plotting
figure(1)
plot(1:length(l1_prob),l1_prob,'-*','Color',color,'LineStyle','-','Marker',symbol)
plot(1:length(l2_prob),l2_prob,'-*','Color',color,'LineStyle',style,'Marker',symbol)

% Plot annotation
grid('on')
% pbaspect([2.5 1 1])
legend_h = legend('location','SouthEast');
set(gca,'XTick',0:10:T_max)
set(gca,'YTick',0:0.1:1)
xlim([T_min T_max])
ylim([0 1])

end

% The New helper function
% This should replace "ProbMatricesWithRank"
% Works? appears so...
function P = PM(m,n,r,q)

P1=zeros(n+1,n+1);
for i=1:length(P1)
    entry_val=1/(q^(n-(i-1)));
    P1(i,i)=entry_val;
    if i<n+1
        P1(i+1,i)=1-entry_val;
    end 
end

s1=zeros(n+1,1);
s1(1)=1;

val=(P1^m)*s1;
P=val(r+1);

end













% Deprecated
% 
% % Prob. of matrices with dimensions (m,n) of rank (r) over field size (q)
% function PMWR = ProbMatricesWithRank(m,n,r,q)
% 
% % The probability og getting rank 0 is = 0 for m=r=0
% if m==0 && r==0
%     PMWR=1;
%     return
% end
% 
% % If m<r there is no possibility to achieve r rank
% if m<r
%     PMWR=0;
%     return
% end
% 
% % Get first set of gaussian coefficients
% gc=gausscoeffs2(n,r,q);
% 
% % Calculate "sum"
% val=0;
% for k=0:r
%     % This should be the one!
%     
%     val2=q^(m*k+binomcoeffs(r-k,2)-n*m);
%     
%     if val2~=0
%         val1=gausscoeffs2(r,k,q);
%         comb_val=val1*val2;
%     else
%         % val2 was zero, so no need to calculate val1
%         comb_val=0;
%     end
%     
%     val=val+((-1)^(r-k)*comb_val);
%     
% end
% 
% % Return probability of matrix 'm'x'n' with rank 'r'
% PMWR=gc*val;
% 
% assert(isnan(PMWR)==0,['ProbMatricesWithRank returned NaN, params: ' num2str([m n r q])]);
% 
% end
% 
% % Compute gaussian coefficients (See Wolfram Alpha)
% % Version 2: Numerical improvements over gausscoeffs(... )
% function GC = gausscoeffs2(m,r,q)
% if r==0
%     % disp('r = 0 in gauss coeffs')
%     GC=1;
% elseif r>0
%     % disp('r > 0 in gauss coeffs')
%     
%     % Calculate numerator
%     % There is X factors in numerator
%     
%     
%     num_length=abs((m))-abs((m-r+1));
%     num=ones(num_length,1);
%     num_index=1;
%     
%     for w=m:-1:m-r+1
%         num(num_index,1)=(q^w-1);
%         num_index=num_index+1;
%     end
%     
%     % Calculate denominator
%     %     denom=1;
%     
%     denom_length=r;
%     denom=ones(denom_length,1);
%     denom_index=1;
%     
%     for w=r:-1:1
%         denom(denom_index,1)=(q^w-1);
%         denom_index=denom_index+1;
%     end
%     
%     % Calculate gaussian coefficient
%     GC=prod(num(:)./denom(:));
%     
%     assert(isnan(GC)==0,'Exception: gausscoeffs2 returned NaN')
%     
% elseif r<0
%     disp('r < 0 error in gausscoeffs!!!')
% end
% 
% end
% 
% % Binomial polynomial thing
% % As on page 123 in "A course in combinatorics"
% function BC = binomcoeffs(a,k)
% 
% tmp_vector=ones(2,1);
% tmp_index=1;
% 
% for w=0:-1:-k+1
%     tmp_vector(tmp_index)=(a+w);
%     tmp_index=tmp_index+1;
% end
% 
% num=prod(tmp_vector);
% denom=factorial(k);
% 
% BC=num/denom;
% 
% end
% 
% 
